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' 1 INTRODUCTION 

v : 

, One of the main longstanding problems of accretion discs in a close binary system is a mechanism of angular momentum 
' transport. Among proposed mechanisms, a-model (Shakura & Sunyaev 1973) is thought to be the standard one. In this model, 

viscosity originating from turbulence, magnetism or whatever, is supposed to transport angular momentum. Nevertheless the 

origin of the viscosity has been poorly understood so far. 

One of other models is the spiral shock model which was first proposed by one of the present authors (Sawada, Matsuda 

6 Hachisu 1986a, b; Sawada et al. 1987). In this model the tidal force due to a companion star excites stationary spiral shock 
waves in the accretion disc. The gas in the disc loses its angular momentum when it passes through the spiral shocks. 

Sawada et al. (1986a, b, 1987) carried out 2D inviscid calculations of a Roche lobe overflow in a semi-detached binary 
system with mass ratio of unity. Their computational region covered both a mass-losing star filling its critical Roche lobe and 
a mass- accreting compact star. They used the second-order-accurate Osher scheme with a generalized curvilinear coordinate. 
They found spiral shocks in an accretion disc for the first time. 

Since then, a number of two-dimensional simulations of accretion discs have been performed mainly with finite differ- 
ence/volume method, and they have confirmed that spiral shocks appear in inviscid discs (Spruit et al. 1987; Spruit 1989; 
Matsuda et al. 1990; Rozyczka & Spruit 1989; Savonije, Papaloizou & Lin 1994). Spruit (1987) obtained self-similar solutions 
of spiral shocks by a semi-analytic manner (see also Chakrabarti 1990). He pointed out that the smaller specific heat ratio 

7 is, the more tightly the spiral arms winds. Godon (1997) studied tidal effects in accretion discs using a two-dimensional 
time-dependent hybrid Fourier- Chebyshev spectral method, under the assumptions of a polytropic equation of state and a 
standard alpha viscosity prescription. He found a critical value of the viscosity (a « 0.01), below which the two-armed spiral 
shock propagates all the way to the inner boundary in cold discs (Mach number M r; 40) for the first time. 

On the other hand, in the case of 3D, this spiral shock model has been criticized on the grounds that the tidally-induced 



ABSTRACT 

We perform 2D and 3D numerical simulations of an accretion disc in a close binary 
system using the Simplified Flux vector Splitting (SFS) finite volume method. In 
our calculations, gas is assumed to be the ideal one, and we calculate the cases with 
7 =1.01, 1.05, 1.1 and 1.2. The mass ratio of the mass-losing star to the mass-accreting 
star is one. Our results show that spiral shocks are formed on the accretion disc in 
all cases. In 2D calculations we find that the smaller 7 is, the more tightly the spiral 
winds. We observe this trend in 3D calculations as well in somewhat weaker sense. 

Recently, Steeghs, Harlaftis & Home (1997) found the first convincing evidence for 
spiral structure in the accretion disc of the eclipsing dwarf nova binary IP Pegasi using 
the technique known as Doppler tomography. We may claim that spiral structure which 
we have discovered in numerical simulations before are now found observationally. 
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perturbations may find it hard to propagate all the way inwards, due to refraction and diversion into the vertical (stratified) 
direction (Lin, Papaloizou & Savonije 1990a, b; Lubow & Pringle 1993). 

Three-dimensional simulations have been carried out mainly using particle methods (Molteni, Belvedere & Lanzafame 
1991; Hirose, Osaki & Mineshige 1991; Nagasawa, Matsuda & Kuwahara 1991; Lanzafame, Belvedere & Molteni 1992, 1993), 
but their results have not shown the spiral shocks. Molteni et al. (1991) and Lanzafame et al. (1992, 1993) stressed that the 
accretion disc itself was not formed for 7 > 1.1 — 1.2. 

Recently, Yukawa, Boffin & Matsuda (1997) made 3D numerical simulations using Smoothed Particle Hydrodynamics 
method (SPH) for three values of 7, i.e. 7=1.01, 1.1 and 1.2, with mass ratio of unity. They used (5 — 6) x 10 4 particles, which 
were much more than those used by previous authors. They used a variable smoothing length technique. They demonstrated 
that the accretion disc as well as spiral shocks existed in the case of 7 = 1.2, contrary to previous claims, although they could 
not obtain clear spiral shocks in the cases of 7=1.1 and 1.01. 

Lanzafame & Belvedere (1997, 1998) performed 3D SPH calculations of an accretion disc in a close binary system via 
wind accretion mechanism. Their simulations were made for a Her X-l-HZ Her like system with a primary mass Mi = 1.3 M 
and a secondary mass M 2 = 2.2 M (1997) and for a Cen X-3 like system with Mi = 1.4 M and M 2 = 19.1 M (1998). 
They considered quasi-polytropic gas (p oc p 1 ) with 7 = 1.01. In their calculation, they calculated the region including a 
portion of the secondary Roche lobe surface. They used 55,747 (1997) and 113,676 (1998) SPH particles. They found that the 
azimuthal distribution of the radial Mach number showed spiral shocks. In a Her X-l-HZ Her case (1997), shocks appeared 
only in the outer regions. On the other hand in a Cen X-3 case (1998), three spiral shocks persisted from the outer edge to 
the inner regions. 

In the case of finite difference/volume method, there had been only one three-dimensional calculation (Sawada & Matsuda 
1992) for a while. They used a Roe upwind TVD scheme to calculate the case of 7 = 1.2 with mass ratio of unity. They 
found the existence of spiral shocks on the accretion disc, but their calculation was done up to only half a revolution period. 
It might be argued that their computational time is not long enough to obtain a convincing result. 

In a series of papers Bisikalo et al. have done 3D numerical studies of the flow structure in semi-detached binary systems 
(1997a, b, 1998a, b, c). They used a TVD Roe scheme modified by monotonic flux limiters in the Osher's form. Their 
computational region was a parallelepipedon [—a. .2a] x [—a. .a] x [0..a], where a is the orbital separation. Non-uniform difference 
grids (finer near the accretor) containing 78 x 60 x 35 grid points for the system X1822-371 and 84 x 65 x 33 grid points for the 
system Z Cha were used (1998b). Their calculations had been stopped at 12-20 revolution periods. In their calculations, gas 
stream from LI did not cause the shock perturbation of the disc boundary. It meant the absence of a hot spot. Their results 
also showed that the gas of 'circumbinary envelope' interacted with the stream and caused the formation of an extended shock 
wave, located on the stream edge. However, their shock waves did not form clear spiral structure. They also claimed that an 
accretion disc was not formed for the case of 7 = 1.2 (1998c) and it was very elongated for 7 = 1.01. As is stated above, not 
only the existence of spiral shocks but that of the disc itself (for higher 7) in 3D has not been established convincingly. 

From observational view points, recently, Steeghs et al. (1997) found the first convincing evidence for spiral structure in 
the accretion disc of the eclipsing dwarf nova binary IP Peg at the outburst phase using the technique known as Doppler 
tomography. 

Neustroev & Borisov (1998) studied the dwarf nova U Geminorum and claimed that they obtained the convincing evidence 
for the existence of spiral shocks in the accretion disc of a cataclysmic variable in quiescence. 

The purpose of the present paper is to perform 2D and 3D finite volume numerical calculations with higher resolution 
up to enough time to confirm the existence of an accretion disc and spiral shocks on it. We also stress that the spirals are 
induced by the tidal force due to the companion. 



2 METHOD OF CALCULATIONS 
2.1 Basic assumptions 

We consider a binary system composed of a mass-accreting primary star with mass Mi and a mass-losing secondary star with 
M2. The mass ratio is defined as q = M2/M1. In the present calculations the mass ratio is fixed to be one. 

We solve only the region around the mass- accreting star. The gas is assumed to be supplied at the LI. In 2D calculations, 
we assume that the flow is confined within a narrow gap in the orbital plane with a constant thickness, so that these calculations 
become two-dimensional. We assume an ideal gas with constant specific heat ratio 7. Although we do not include heating and 
cooling effects explicitly, we choose the smaller 7 so that these effects are implicitly taken into account. In the present paper, 
we examine the cases with 7=1.2, 1.1, 1.05 and 1.01. These calculations contain no viscosity except numerical dissipation. 

In our 3D calculations, the assumptions of calculations are the same as in 2D calculations except that gas is not restricted 
within a gap of constant thickness. Symmetry about the orbital plane is assumed. 
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2.2 Basic equations 



The basic equations we are going to solve are the two-dimensional as well as the three-dimensional Euler equations describing 
the motion of an inviscid gas. 

We choose the orbital separation, a, as the length scale, and a uj as the velocity scale, where U) is the angular velocity of 
the system. Therefore, the revolution period is normalized to 2n. Density is normalized by the density of the gas at LI point. 
The primary star is centred at the origin and the secondary star at (-1, 0, 0). We set a numerical grid centred at the origin 
and just touching the LI point. Cartesian coordinates we use have 200 x 200 grid points in 2D and 200 x 200 x 50 grid points 
in 3D, respectively. 

Our 3D basic equations in the rotating frame of reference can be described in dimensionless form as 
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where p,p, e, u, v, w, k x , k v , and k z are the density, the pressure, the total energy per unit volume, the x, y and z components 
of velocity and the x, y, and z components of force, respectively. The x, y and z components of force are 
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Here, rt(i — 1, 2) is the distance from the point considered to the centre of each star, and Xi, yi and Zi are each star's Cartesian 
components, ig represents the x component of the distance from the centre of mass. 

In our calculations, the mass-accreting star is centred at the origin, and the mass-losing star at (-1, 0, 0). Hence 



Xi = x, 
X2 = X + 1 
XG = X + 
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1 + 9 

This geometry is shown in Fig. |l[ 

In equations and (Q), the first term represents the Coriolis force, the second one the centrifugal force, respectively. 
The last two terms in equations (^) and (^]) as well as ones in equation Q represent the gravity force. 

The equation of state we use is 



Pi 2 . 2 2n 
e — —(u +v + w ) 



P= (7- 1) 

In the normalization procedure, we use an auxiliary equation: 
aV=G(Mi+M 2 ), 

to eliminate the gravitational constant G from the equations. 



(9) 



(10) 



2.3 Numerical method 

We use Simplified Flux vector Splitting (SFS) scheme (Jyounouti et al. 1993; Shima & Jyounouti 1994, 1997) (see Appendix). 
In order to obtain high spatial resolution we use a MUSCL (Monotone Upwind Scheme for Conservation Law) type approach. 
The numerical accuracy of the method is the second-order both in space and in time. 



4 Makoto Makita, Kenji Miyawaki and Takuya Matsuda 



y 
A 





o [mT] 



»- x 



Figure 1. Geometry of the system. 



2.4 Initial and boundary conditions 

As an initial condition we suppose the whole computational region, which is —0.5 < x,y < 0.5 in 2D and furthermore < 
z < 0.25 in 3D, is filled with tenuous gas with higher temperature at the initial instance. We assume po = 10~ 5 , po = 10~ 4 /7 
and uo = vo = wq = in the whole numerical region. 

This initial condition is also maintained outside of the outer boundary except for LI point and inside of the inner boundary 
during calculations. This choice forms the outer and the inner boundary conditions except the inlet at LI point. 

At LI point a small rectangular inlet hole is placed. The values at the inlet are p in — 1.0, pt„ = 10~ 2 /7, u in = 0.01 
and v in — Wi n = 0.0 at all the time. This choice of the density and the pressure means that the sound speed of the gas to 
be injected into the computational region is 10 _1 . Therefore, we may say that the temperature of the gas is very high. This 
is necessary to ensure a sufficient amount of gas flux from the inlet. The mass flux from LI point is computed by solving a 
Ricmann problem, so the amount of inflow gas is not the same between 2D and 3D simulations. 

The mass-accreting star is represented by a hole of 3 x 3 cells in 2D cases and 3x3x2 cells in 3D cases situated at the 
origin. We follow the time evolution until t = 207r (that is 10 revolution periods) in 2D calculations and t = 10ir in 3D cases. 

3 RESULTS 

3.1 2D calculations 

The time evolution of the disc mass is shown in Fig. ^. In the case of 7 = 1.2, the mass is saturated at about t = 50. The disc 
mass is still increasing except for the case of 7 = 1.2. We have not yet reached a steady state in the other cases. In order to 
investigate a time variation of flow, we make video movies. From the movies we confirm that a gross feature of the flow in our 
calculations does not change except minor oscillations. We are confident that we have almost obtained the final structure of 
the disc. 

We show the gray scale of density distribution for all cases at t = 44, which is 7 revolution periods, in Fig. [| We find 
two clear spiral shocks in these figures as the same as the results of the previous authors. The inflow from LI point expands 
rapidly because of rather high temperature of the gas assumed. This flow is a kind of an under-expanded jet; we observe two 
intercepting shocks are formed. These shocks are artifact due to our boundary condition at LI point. In order to remove these 
shocks, we have to solve the whole region containing the mass-losing star as was done by Sawada et al. (1986a, b, 1987). 
However, we believe that the structure of the accretion disc is not affected by the present choice of the boundary condition 
at LI point. It can be proved by comparing the present results with our previous ones (Sawada et al. 1986a, 1987). 

From the point of view of dependence on 7, we find that the spiral arms wind more tightly in smaller 7 cases than larger 
ones. Lower 7 means a cooler disc with larger Mach number of the flow. Our results coincide with those of Savonije et al. 
(1994) and Godon (1997), and the analytical work of Spruit (1987). Spruit also argued that eventually spiral waves become 
a ring for 7=1. Our results show that spiral shocks still exist even for 7 = 1.01. However, our disc has a constant thickness, 
while Spruit considered a tapered disc. Therefore, direct comparison may not be appropriate. 

We also observe that the flow for higher 7 are more unsteady in time, of which fact may be observed in the case of 7 = 1.2 
in Fig. H and also in our video movies. 
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Figure 2. Evolution of the disc mass as a function of time. 

3.2 3D calculation 

We present density profiles in the orbital plane in all four cases in Fig. ^. These figures show the existence of spiral shocks 
in not only the case of 7 = 1.2 but also in the other three cases. We note the remarkable similarity between the present 
case of 7 = 1.2 and those by Sawada et al. (1992) and by Yukawa et al. (1997). However, Yukawa et al. (1997) could not 
obtain spiral structures for 7 = 1.1 and 1.01. We may attribute this to the lack of number of particles in their calculations. 
As to the accuracy of the computational results, we would point out that our number of grid points is 4 x 10 6 for the whole 
computational domain compared with (5 — 6) x 10 4 SPH particles in Yukawa et al. (1997). Although it may not be appropriate 
to compare these numbers directly, the difference of two orders of magnitude in these numbers may be significant. In fact 
Lanzafame & Belvedere (1998) observed spiral structures for the case of 7 = 1.01 using their SPH code with 113,676 particles. 
Therefore, we believe that there should exist spiral structures on the accretion discs as long as 7 < 1.2. 

Comparing our case of 7 = 1.01 with our case of 7 = 1.2, we see that the former case shows slightly tightly wound spiral 
shocks near the centre. We, however, do not find a clear difference in spiral structure in the outer region particularly for the 
cases of 7 < 1.1. We, therefore, may conclude that the dependence of spiral structure on the value of 7 in 3D calculation is 
somewhat weaker than that in the 2D cases. 

We make video movies for all cases and confirm that the flow reaches almost steady state. It is observed that the higher 
7 cases are more violent than the lower 7 cases, of which fact is the same as in 2D cases. 

We provide the three-dimensional views of the iso-density surface of p = 10~ 3 in Fig. ^. Fig. [| may lead an impression 
that the stream from LI point goes into the disc directly without the interaction between the stream and the disc. However, 
this is certainly not the case. Fig. ^ clearly shows that the stream does not penetrate into but collides with the disc. Therefore 
we do not confirm the claim by Bisikalo et al. (1998b, c) on the non-existence of a so-called hot spot. 

In Fig. |B| we observe the accretion disc, the spiral structures in it and the stream from LI point. We stress here again 
that the accretion disc does exist for the case of 7 = 1.2 contrary to the claims by Lanzafame et al. (1992) and Bisikalo et al. 
(1998c). Our result is consistent with our earlier finite volume calculation by Sawada et al. (1992) and the SPH calculation 
by Yukawa et al. (1997). 

The strange shape of the stream from LI point is due to the fact that the shape of the inlet at the LI point is rectangular 
rather than circular, and the stream is an under expanded jet. Nevertheless, it is important to stress again that the stream is 
an independent object from the accretion disc. 

The cause of the discrepancy between our result and Bisikalo et al. (1998b, c) is not well understood. However, we should 
point out that our number of grid points above the orbital plane is 2 x 10 6 , while theirs is at most 1.8 x 10 5 , i.e. ours is more 
than 10 times more than theirs. Moreover the volume of our computational domain is 0.25, while theirs is 6. It means that 
the volume of their computational cell is almost 240 times bigger than ours on average, which suggests that their calculations 
are much cruder than ours. We are not sure if this is the sole reason of the discrepancy, though. 
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Figure 3. Gray scale of density distribution with logarithmic scale at 7 revolution periods. Bar in the right side shows the scale range. 

Fig. |i| presents the density profiles in the x-z plane. We can see that disc thickness near the companion star is thicker 
than the other side. This feature coincides with the result of Yukawa et al. (1997). We also find that the smaller 7 is, the 
more symmetry the disc structure is. 

We make a test calculation in which gas inflow is terminated at t = 10-7T to see the effect of the inflow on the spiral shocks. 
We show the density profiles in the case of 7 = 1.1 at one revolution period after the inflow is terminated, in Fig. [?]. We 
observe clear two-armed spiral shocks of bi-symmetry. The spiral shocks are also found even after two revolution period. We, 
therefore, conclude that the spiral shocks is formed by tidal force, not by the inflow, of which claim was posed by Bisikalo et 
al. (1998b, c). 

Figs. ^ and [] show Mach number contours with gray scale as well as the velocity vectors in the orbital plane at 5 revolution 
periods. In these figures, the velocity vectors change the direction abruptly at the positions corresponding to the spirals shown 
in Fig. ^. We find also that Mach number in accretion discs become higher for smaller 7. This feature coincides with the 
results of previous 2D calculations. It is clear that the velocity in the outer region is nearly circular and slower for all cases. 
This feature in 7 = 1.2 case also contradicts very much with Bisikalo et al. (1998c), who found a large mass outflow from the 
system. 

To confirm the existence of spiral shocks and to compare with the observations by Steeghs et al. (1997), we make a 
Doppler map, which is essentially a hodograph used in hydrodynamics. Doppler map shows the distribution of a physical 
variable in velocity coordinates (V x , V v ). The relation between the Doppler map and the physical space is as follows: 
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Figure 4. Gray scale of density distribution with logarithmic scale in the orbital plane at t = 27. 



Figure 5. 3D view of iso-density surface at the level of p = —3 with logarithmic scale at 5 revolution periods. Gray scale of density in 
the orbital plane is also shown. Bar in the left side shows the scale range. 



(i) Doppler map is 90° rotated anti-clockwise (in our case) from the physical space. 

(ii) The inner region of the accretion disc, which has faster rotational velocity, is mapped to the outer region of the Doppler 
map, and vice versa. 

(iii) Spiral structure in the physical space is mapped to that in the Doppler map. 

From the property (iii) we may say that, if spiral structure is observed in a Doppler map, it shows the existence of spiral 
structure in a physical space. (See Marsh & Home (1988), Steeghs et al. (1996) and Steeghs et al. (1997)). This is the basis 
of their claim of the discovery of spiral structure in an accretion disc in IP Peg by Steeghs et al. (1997). 

Figs. and show density distribution in Doppler coordinates (V x , V y ) for 2D and 3D cases. From these maps, we find 
the clear spiral structure in both 2D and 3D (and the stream from LI point.) We, therefore, conclude that there exists the 
spiral structure in real physical space. Since the mass ratio of the binary in IP Peg is not 1 but 0.5, the direct comparison 
between Figs. 10 or 11 with the observation may be meaningless. We calculated 2D flows for the mass ratio of 0.5 and compared 
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Figure 6. Gray scale of density distribution in the x-z plane at t = 27. The scale range is logarithmic one which is the same as Fig. 



with the observation (Matsuda et al. 1998) and obtained a reasonable agreement. The comparison between 3D calculations 
and the observations is now under way. 



4 CONCLUSION 

We have performed 2D and 3D numerical simulations of accretion discs in close binary systems using SFS finite volume 
method. We obtain the following conclusions: 

(i) In 2D calculations we observe clear spiral shocks, which winds more tightly for smaller 7. (See Fig. This result is 
consistent with previous works. 

(ii) In 3D cases we confirm that discs are formed even for 7 > 1.1 contrary to the claims by Lanzafame et al. (1992) and 
Bisikalo et al. (1998c) (see Fig. ||), but this is consistent with the previous works (Sawada & Matsuda 1992; Yukawa et al. 
1997). 

(iii) In 3D calculations we also find spiral shocks for all four cases, i.e. 7 = 1.2, 1.1, 1.05 and 1.01, as well as 2D cases. 

(iv) However, the pitch angle of the spirals in 3D is not so markedly correlated with 7 as in 2D cases. (See Figs. | and ym. 

(v) The stream from the LI point does not smoothly penetrate into the accretion disc, of which fact contradicts with the 
claim by Bisikalo et al. (1998b, c). (See Fig. |). 

(vi) Spiral Shocks are produced by the tidal interaction contrary to the claim by Bisikalo et al. (1998b, c). (See Fig. Q). 

(vii) Doppler maps for 2D and 3D are constructed. Such maps can be compared with observations. 
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Figure 8. Mach number contours with gray scale and velocity vectors in the orbital plane at 5 revolution periods. Top: The case of 
7 = 1.2, Bottom: The case of 7 = 1.1. 
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APPENDIX A: SFS SCHEME 

SFS scheme is one of Advection Upstream Splitting Method (AUSM) type schemes. Original AUSM scheme invented by Liou 
& Steffen (1993) is very simple, robust for strong shock and accurate, but shows overshoot at shock front. Inspired by AUSM, 
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Figure 9. Mach number contours with gray scale and velocity vectors in the orbital plane at 5 revolution periods. Top: The case of 
7 = 1.05, Bottom: The case of 7 = 1.01. 



various schemes including SFS have been proposed since then. Although, various schemes have been developed independently, 
these schemes can be written in a common form. They can be called as AUSM type schemes. 

Now we show the formulation of AUSM type scheme for three-dimensional Euler equation. (See Shima & Jounouchi 
1997). Three-dimensional Euler equation can be written in the integral form as follows: 



Fds = 0, 



(Al) 



Q = 



( p \ 

pu 
pv 
pw 

V e J 



(A2) 
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Figure 10. Doppler maps (hodograph) of 2D calculations in the orbital plane at 5 revolution periods. Density profiles in velocity 
coordinates (V x ,V y ) are shown. 



( 1 \ 



F = m* + pN, * 
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m = pV n , V n = x n u + y n v + z n w, 



(A3) 



(A4) 



where p, u, v, w, e,p, and h — (e + p)/p represent the density, the velocity in the x, y and z direction, the total energy per unit 
volume, the pressure and the total enthalpy, respectively. x n , y n and z„ show unit normals to the surface. This form suggests 
that Euler flux can be divided into the advection term and the pressure term. AUSM is derived based on the fact that the 
advection term and the pressure term can be upwinded separately. 
AUSM type schemes can be written in the following form, 



p _ M + 



i*_ +pN, 



(A5) 



where subscript ± show physical value at left (+) and right (— ) side of a cell boundary, and p does mixing of pressure using 
Mach number of left and right state which is defined by 



P = P+P+ + P-P- 



(A6) 
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Figure 11. Doppler maps (hodograph) of 3D calculations in the orbital plane at 5 revolution periods. Density profiles in velocity 
coordinates (V x ,V y ) are shown. 

P± = ^(2 T M±)(M±±l) 2 , if \M±\<1, (A7) 

and f}± are smoothly switched to 1 or for supersonic case. 

By choosing different forms for the mass flux m above, various AUSM type schemes have been proposed. SFS uses 
variations of van Leer's Flux Vector Splitting (FVS) for its mass flux. This is written in following form, 

m = m + + m_ , (A8) 

m± = ±£±A( M ±±l) 2 , if\M±\<l, (A9) 
4c 

M± = c± = — , c± = 7 — , (A10) 

c± c p± 

where 



-T- (An) 

c is arithmetic average of sound speed and pure upwind side value or of M± will be used for supersonic case. SFS has 
improved overshoot at a shock. 

This paper has been produced using the Royal Astronomical Society/Blackwell Science I^TpjX style file. 
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